
S = 0.3
C = [0.05,0.1,0.2]

n_A = 51L
A = dindgen(n_A)/100.d

n_Td = 81L
Td = dindgen(n_Td)*0.25d + 10.d
print,Td

ffore = fltarr(n_A,n_Td)
restore,'./irdc_dist_model/data/case3_irdc_grid_AT.sav',/ver

myps,'./irdc_dist_model/analysis_plots/ffore_A_Td.eps',xsize=14,/cmyk,ysize=9
multiplot,[3,2],/square,/doxaxis,$;mytitle='Dust Temperature [K]',$
          mxtitle='Foreground "veil" (A)',ygap=0.04,mtitoff=2,$
          mtitle='Hypothetical BGPS Source with S!d40!n = '+string(S,format="(F0.2)")+' Jy'

;; Loop over Contrasts
FOR i=0,2 DO BEGIN
   print,'Next Plot'
   ;; Build Ffore
   FOR j=0,n_A-1 DO BEGIN
      FOR k=0,n_Td-1 DO BEGIN
         
         ffore[j,k] = -A[j] + 1. - (C[i] / (1. - exp(-flux2tau_bgps(Td[k])*S)))
         
      ENDFOR
   ENDFOR
   
   xr = (i EQ 2) ? [0,.5] : [0,.499999]
   ytf = (i EQ 0) ? "(I0)" : 'blank_axis'
   ytit = (i EQ 0) ? 'Dust Temperature [K]' : ''
   print,ytf
   
   loadct,13,/silent
   
   print,m4_Stat(ffore)
   plotimage,ffore,range=[0,1],xst=4,yst=4,ytit=ytit,$
             title='Contrast = '+string(C[i],format="(F0.2)")
   axis,xaxis=0,xr=xr,/xst,/save
   axis,xaxis=1,xr=xr,/xst,/save,xtickformat='blank_axis'
   axis,yaxis=0,yr=[10,30],/yst,/save,ytickformat=ytf
   axis,yaxis=1,yr=[10,30],/yst,/save,ytickformat='blank_axis'
   
   levels = min(grid) + [2.3, 6.2, 11.8, 20, 30, 40, 50, 100]
   colors = [replicate('WT3',3),replicate('RED6',4),'WT3']
   print,levels
   gTd = dindgen(251)*0.1d + 10.d
   cgContour,grid,A,gTd,levels=levels,/overplot,c_colors=colors,$
           c_thick=5,c_labels=0
   
   ;; cgContour,ffore,A,Td,levels=[0,0.2,0.4,0.6,0.8,1.0],c_thick=3,/over,$
   ;;           c_colors=[replicate('White',4),replicate('Black',2)],$
   ;;           c_charsize=0.8
   
   loadct,0,/silent
   multiplot,/doyaxis,/doxaxis
ENDFOR

loadct,13,/silent
cgColorbar,range=[0.,1.],/vertical,/right,position=[0.94,0.57,0.96,0.91],$
           title='f!dfore!n',charsize=1.0
;; /vert,/right,position=[0.92,0.18,0.95,0.90],$
;;            title='f!dfore!n',range=[0,1],charsize=1.0
;; ;,format="(F0.1)",charsize=1.0


;; myps,/done,/mp

;;====================================================================
;; Flip Side

ff = [0.2,0.5,0.8]
contrast = fltarr(n_A,n_Td)

;; myps,'./irdc_dist_model/analysis_plots/contrast_A_Td.eps',xsize=14,/cmyk
;; multiplot,[3,1],/square,mytitle='Dust Temperature [K]',$
;;           mxtitle='Foreground "veil" (A)'

;; Loop over Contrasts
FOR i=0,2 DO BEGIN
   print,'Next Plot'
   ;; Build Contrast
   FOR j=0,n_A-1 DO BEGIN
      FOR k=0,n_Td-1 DO BEGIN
         
         contrast[j,k] = (-A[j]+1.-ff[i])*(1. - exp(-flux2tau_bgps(Td[k])*S))
         
      ENDFOR
   ENDFOR
   
   xr = (i EQ 2) ? [0,.5] : [0,.499999]
   ytf = (i EQ 0) ? "(I0)" : 'blank_axis'
   ytit = (i EQ 0) ? 'Dust Temperature [K]' : ''
   print,ytf
   
   loadct,13,/silent
   
   print,m4_Stat(ffore)
   pr = [0,0.3]
   plotimage,contrast,range=pr,xst=4,yst=4,ytit=ytit,$
             title='f!dfore!n = '+string(ff[i],format="(F0.2)")
   axis,xaxis=0,xr=xr,/xst,/save
   axis,xaxis=1,xr=xr,/xst,/save,xtickformat='blank_axis'
   axis,yaxis=0,yr=[10,30],/yst,/save,ytickformat=ytf
   axis,yaxis=1,yr=[10,30],/yst,/save,ytickformat='blank_axis'
   
   levels = min(grid) + [2.3, 6.2, 11.8, 20, 30, 40, 50, 100]
   colors = [replicate('WT3',3),replicate('RED6',4),'WT3']
   print,levels
   gTd = dindgen(251)*0.1d + 10.d
   cgContour,grid,A,gTd,levels=levels,/overplot,c_colors=colors,$
           c_thick=5,c_labels=0
   
   ;; cgContour,ffore,A,Td,levels=[0,0.2,0.4,0.6,0.8,1.0],c_thick=3,/over,$
   ;;           c_colors=[replicate('White',4),replicate('Black',2)],$
   ;;           c_charsize=0.8
   
   loadct,0,/silent
   multiplot,/doyaxis
ENDFOR

loadct,13,/silent
cgColorbar,range=pr,/vertical,/right,position=[0.94,0.14,0.96,0.48],$
           title='Contrast',charsize=1.0
;; /vert,/right,position=[0.92,0.18,0.95,0.90],$
;;            title='f!dfore!n',range=[0,1],charsize=1.0
;; ;,format="(F0.1)",charsize=1.0


myps,/done,/mp
END
